model <- ModelTBBCGEngland::read_libbi(file.path(model_dir, "libbi", "posterior"))
priors <- readRDS(file.path(model_dir, "data", "prior-params.rds"))

model
## Wrapper around LibBi
## ======================
## Model:  Baseline 
## Run time:  3790.562  seconds
## Number of samples:  1000 
## State trajectories recorded:  E foi H L N P S T_E T_P YearlyAgeCases YearlyCases YearlyEAgeCases YearlyECases YearlyEPulCases YearlyEPulDeaths YearlyNonUKborn YearlyPAgeCases YearlyPCases YearlyPulCases YearlyPulDeaths 
## Observation trajectories recorded:  YearlyAgeInc YearlyHistPInc YearlyInc 
## Parameters recorded:  alpha_t c_eff c_hist chi_init delta epsilon_h_0_4 epsilon_h_15_89 epsilon_h_5_14 epsilon_l_0_4 epsilon_l_15_89 epsilon_l_5_14 HistMeasError kappa_0_4 kappa_15_89 kappa_5_14 M mu_e_0_14 mu_e_15_59 mu_e_60_89 mu_p_0_14 mu_p_15_59 mu_p_60_89 nu_e_0_14 nu_e_15_89 nu_p_0_14 nu_p_15_89 phi_0_14 phi_15_59 phi_60_89 rho_0_14 rho_15_59 rho_60_89 Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89 zeta_p_0_14 zeta_p_15_59 zeta_p_60_89
traces <- coda::mcmc(rbi::get_traces(model))

coda::rejectionRate(traces)
##               M           c_eff          c_hist           delta 
##        0.989990        0.989990        0.989990        0.989990 
##   epsilon_h_0_4  epsilon_h_5_14 epsilon_h_15_89       kappa_0_4 
##        0.989990        0.989990        0.989990        0.989990 
##      kappa_5_14     kappa_15_89   epsilon_l_0_4  epsilon_l_5_14 
##        0.989990        0.989990        0.989990        0.989990 
## epsilon_l_15_89        phi_0_14       phi_15_59       phi_60_89 
##        0.989990        0.989990        0.989990        0.989990 
##    Upsilon_0_14   Upsilon_15_59   Upsilon_60_89        rho_0_14 
##        0.990991        0.989990        0.990991        0.989990 
##       rho_15_59       rho_60_89       nu_p_0_14      nu_p_15_89 
##        0.989990        0.989990        0.989990        0.989990 
##       nu_e_0_14      nu_e_15_89     zeta_p_0_14    zeta_p_15_59 
##        0.989990        0.989990        0.989990        0.989990 
##    zeta_p_60_89     zeta_e_0_14    zeta_e_15_59    zeta_e_60_89 
##        0.989990        0.989990        0.989990        0.989990 
##       mu_p_0_14      mu_p_15_59      mu_p_60_89       mu_e_0_14 
##        0.989990        0.989990        0.989990        0.989990 
##      mu_e_15_59      mu_e_60_89        chi_init       alpha_t.0 
##        0.989990        0.989990        0.989990        1.000000 
##       alpha_t.1       alpha_t.2       alpha_t.3       alpha_t.4 
##        1.000000        1.000000        1.000000        1.000000 
##       alpha_t.5   HistMeasError 
##        1.000000        0.989990
plot_param(model, scales = "free", plot_type = "trace") + theme(legend.position = "none")
## Aggregate function missing, defaulting to 'length'

- Plot overview of prior and posterior densities

plot_param(model, prior_params = priors, scales = "free")
## Aggregate function missing, defaulting to 'length'

param_sum <- plot_param(model, prior_params = priors, plot_data = FALSE) %>% 
  group_by(Distribution, parameter, length) %>% 
  summarise(median = median(value), 
            lll = quantile(value, 0.025), 
            hhh = quantile(value, 0.975)) %>% 
  mutate(value = pretty_ci(median, lll, hhh, sep = ", ")) %>%
  group_by(Distribution, parameter) %>% 
  mutate(Parameter = case_when(max(length) > 1 ~ paste(parameter, 1:n(), sep = "-"),
                               TRUE ~ as.character(parameter))
         ) %>% 
  ungroup %>% 
  select(Distribution, Parameter, value) %>% 
  spread(key = "Distribution", value = "value") %>% 
  select(Parameter, Prior, Posterior)
## Aggregate function missing, defaulting to 'length'
saveRDS(param_sum, file.path(model_dir, "data", "sum_prior_posterior.rds"))
knitr::kable(param_sum)
Parameter Prior Posterior
alpha_t-1 0.85 (0.77, 0.90) 0.86 (0.86, 0.86)
alpha_t-2 0.69 (0.52, 0.80) 0.75 (0.75, 0.75)
alpha_t-3 0.57 (0.33, 0.71) 0.45 (0.45, 0.45)
alpha_t-4 0.57 (0.37, 0.71) 0.54 (0.54, 0.54)
alpha_t-5 0.25 (-0.11, 0.48) 0.09 (0.09, 0.09)
alpha_t-6 0.21 (-0.44, 0.55) 0.14 (0.14, 0.14)
c_eff 2.61 (0.14, 4.85) 4.87 (4.87, 4.87)
c_hist 12.48 (10.16, 14.87) 10.82 (10.82, 10.82)
chi_init 0.19 (0.08, 0.28) 0.25 (0.25, 0.25)
delta 0.78 (0.70, 0.87) 0.70 (0.70, 0.70)
epsilon_h_0_4 1.49 (0.09, 5.16) 4.40 (4.40, 4.40)
epsilon_h_15_89 24.65 (1.42, 73.48) 94.34 (94.34, 94.38)
epsilon_h_5_14 3.90 (0.18, 11.83) 9.50 (9.50, 9.51)
epsilon_l_0_4 616.42 (33.73, 1720.16) 51.48 (51.41, 51.52)
epsilon_l_15_89 1081.98 (58.39, 3378.35) 1766.49 (1765.83, 1766.49)
epsilon_l_5_14 489.66 (28.10, 1488.47) 921.11 (921.11, 922.35)
HistMeasError 1.00 (0.61, 1.40) 1.25 (1.25, 1.25)
kappa_0_4 0.86 (0.05, 2.79) 0.67 (0.67, 0.67)
kappa_15_89 1.09 (0.05, 3.41) 2.90 (2.90, 2.90)
kappa_5_14 0.97 (0.06, 3.10) 1.44 (1.44, 1.44)
M 0.24 (0.01, 0.49) 0.49 (0.49, 0.49)
mu_e_0_14 275.68 (217.01, 341.37) 239.52 (239.52, 239.52)
mu_e_15_59 193.36 (72.76, 318.29) 401.85 (401.72, 401.85)
mu_e_60_89 37.15 (2.08, 103.08) 71.24 (71.24, 71.31)
mu_p_0_14 241.67 (155.25, 324.07) 190.08 (190.07, 190.08)
mu_p_15_59 82.93 (4.33, 264.25) 248.62 (248.56, 248.62)
mu_p_60_89 50.54 (3.38, 158.42) 39.28 (39.27, 39.28)
nu_e_0_14 0.19 (0.00, 1.16) 0.19 (0.19, 0.19)
nu_e_15_89 0.31 (0.01, 1.93) 1.15 (1.15, 1.15)
nu_p_0_14 0.11 (0.00, 0.78) 0.01 (0.01, 0.01)
nu_p_15_89 0.24 (0.01, 1.25) 4.40 (4.40, 4.40)
phi_0_14 0.57 (0.28, 1.04) 0.53 (0.53, 0.53)
phi_15_59 0.61 (0.26, 1.18) 1.72 (1.72, 1.72)
phi_60_89 0.59 (0.26, 1.11) 0.50 (0.50, 0.50)
rho_0_14 0.30 (0.26, 0.34) 0.27 (0.27, 0.27)
rho_15_59 0.65 (0.64, 0.66) 0.69 (0.69, 0.69)
rho_60_89 0.54 (0.52, 0.55) 0.56 (0.56, 0.56)
Upsilon_0_14 0.63 (0.61, 0.65) 0.60 (0.60, 0.60)
Upsilon_15_59 0.71 (0.70, 0.71) 0.70 (0.70, 0.70)
Upsilon_60_89 0.75 (0.74, 0.76) 0.76 (0.76, 0.76)
zeta_e_0_14 122.39 (57.91, 187.56) 142.61 (142.56, 142.61)
zeta_e_15_59 60.81 (3.67, 171.01) 41.51 (41.51, 41.58)
zeta_e_60_89 135.73 (55.90, 212.67) 98.79 (98.78, 98.85)
zeta_p_0_14 94.95 (18.70, 179.68) 19.81 (19.81, 19.89)
zeta_p_15_59 78.52 (3.22, 242.86) 605.28 (605.23, 605.35)
zeta_p_60_89 121.96 (15.39, 248.56) 63.16 (63.13, 63.16)
plot_state(model, summarise = TRUE)

plot_state(model, summarise = TRUE, start_time = 59)

plot_state(model, states = c("YearlyHistPInc", "YearlyInc"), summarise = TRUE, start_time = 59)

pred_obs <- plot_state(model, states = c("YearlyHistPInc", "YearlyInc", "YearlyAgeInc"), summarise = FALSE, start_time = 59, plot_data = FALSE)

obs <- ModelTBBCGEngland::setup_model_obs() %>% 
  bind_rows(.id = "state")

pred_obs <- pred_obs %>% 
  dplyr::filter(Average == "median") %>% 
  mutate(pred_value = pretty_ci(Count, lll, hhh, digits = 0)) %>% 
  left_join(obs, by = c("state" = "state", "time" = "time", "age" = "age")) %>% 
  select(Observation = state, Age = age, Year = time,  `Observed Incidence` = value, `Predicted Incidence` = pred_value) %>% 
  mutate(Year = Year + 1931)


sum_obs_table <- pred_obs %>% 
  dplyr::filter(Observation %in% c("YearlyInc", "YearlyHistPInc")) %>% 
  drop_na(`Observed Incidence`) %>% 
  select(-Age) %>% 
  mutate(Observation = case_when(Observation == "YearlyInc" ~ "All TB cases",
                                 Observation == "YearlyHistPInc" ~ "Pulmonary TB cases"))


saveRDS(sum_obs_table, file.path(model_dir, "data", "sum_obs_table.rds"))

age_cases_table <- pred_obs %>% 
  dplyr::filter(Observation %in% c("YearlyAgeInc")) %>% 
  drop_na(`Observed Incidence`) %>% 
  select(-Observation) %>% 
  group_by(Year) %>% 
  mutate(Age = c(paste0(seq(0, 65, 5), "-", seq(4, 69, 5)), "70-89")) %>% 
  ungroup


saveRDS(age_cases_table, file.path(model_dir, "data", "age_cases_table.rds"))
knitr::kable(sum_obs_table)
Observation Year Observed Incidence Predicted Incidence
Pulmonary TB cases 1990 3618 4154 (4147 to 4366)
Pulmonary TB cases 1991 3596 4063 (4062 to 4106)
Pulmonary TB cases 1992 3816 3804 (3804 to 4009)
Pulmonary TB cases 1993 3961 3748 (3748 to 3877)
Pulmonary TB cases 1994 3694 3657 (3608 to 3806)
Pulmonary TB cases 1995 3711 3608 (3595 to 3723)
Pulmonary TB cases 1996 3690 3581 (3544 to 3700)
Pulmonary TB cases 1997 3947 3632 (3566 to 3717)
Pulmonary TB cases 1998 3926 3466 (3466 to 3594)
Pulmonary TB cases 1999 3827 3506 (3476 to 3578)
All TB cases 2000 1803 1811 (1774 to 1937)
All TB cases 2001 1866 1817 (1766 to 1817)
All TB cases 2002 1833 1779 (1684 to 1784)
All TB cases 2003 1685 1652 (1612 to 1772)
All TB cases 2004 1776 1688 (1669 to 1749)
knitr::kable(age_cases_table)
Age Year Observed Incidence Predicted Incidence
0-4 2000 75 95 (79 to 98)
5-9 2000 59 55 (55 to 74)
10-14 2000 75 74 (64 to 77)
15-19 2000 112 14 (14 to 23)
20-24 2000 133 35 (26 to 36)
25-29 2000 116 60 (54 to 66)
30-34 2000 147 75 (60 to 90)
35-39 2000 99 121 (102 to 121)
40-44 2000 83 146 (143 to 166)
45-49 2000 105 192 (155 to 192)
50-54 2000 116 204 (189 to 204)
55-59 2000 112 228 (220 to 228)
60-64 2000 101 228 (209 to 236)
65-69 2000 99 182 (182 to 214)
70-89 2000 371 137 (137 to 145)
0-4 2001 113 84 (78 to 84)
5-9 2001 65 52 (52 to 67)
10-14 2001 51 82 (70 to 89)
15-19 2001 130 19 (13 to 24)
20-24 2001 145 41 (34 to 44)
25-29 2001 127 46 (46 to 72)
30-34 2001 122 71 (71 to 92)
35-39 2001 128 111 (96 to 127)
40-44 2001 107 152 (128 to 164)
45-49 2001 101 185 (175 to 194)
50-54 2001 96 213 (198 to 221)
55-59 2001 91 232 (206 to 233)
60-64 2001 94 214 (184 to 224)
65-69 2001 105 198 (196 to 235)
70-89 2001 391 151 (121 to 151)
0-4 2002 102 85 (69 to 100)
5-9 2002 62 67 (54 to 71)
10-14 2002 64 75 (66 to 75)
15-19 2002 129 16 (12 to 17)
20-24 2002 148 37 (25 to 37)
25-29 2002 120 46 (44 to 55)
30-34 2002 112 76 (64 to 83)
35-39 2002 127 86 (86 to 115)
40-44 2002 98 138 (127 to 151)
45-49 2002 89 175 (159 to 176)
50-54 2002 104 189 (189 to 219)
55-59 2002 116 203 (203 to 217)
60-64 2002 88 215 (200 to 223)
65-69 2002 90 195 (190 to 200)
70-89 2002 384 145 (131 to 159)
0-4 2003 74 102 (80 to 102)
5-9 2003 46 51 (51 to 71)
10-14 2003 59 55 (55 to 72)
15-19 2003 102 15 (11 to 17)
20-24 2003 116 24 (24 to 28)
25-29 2003 137 54 (49 to 54)
30-34 2003 118 66 (66 to 73)
35-39 2003 113 108 (103 to 108)
40-44 2003 109 115 (115 to 149)
45-49 2003 83 185 (153 to 185)
50-54 2003 102 188 (184 to 200)
55-59 2003 92 203 (190 to 212)
60-64 2003 98 203 (184 to 208)
65-69 2003 94 224 (182 to 224)
70-89 2003 342 126 (126 to 145)
0-4 2004 112 106 (68 to 106)
5-9 2004 71 73 (57 to 76)
10-14 2004 81 98 (72 to 98)
15-19 2004 111 19 (14 to 19)
20-24 2004 120 35 (11 to 35)
25-29 2004 136 37 (37 to 48)
30-34 2004 129 45 (45 to 71)
35-39 2004 103 99 (87 to 99)
40-44 2004 128 120 (110 to 137)
45-49 2004 94 157 (148 to 186)
50-54 2004 103 177 (177 to 196)
55-59 2004 98 221 (195 to 221)
60-64 2004 86 214 (193 to 214)
65-69 2004 92 194 (192 to 202)
70-89 2004 312 150 (131 to 150)
plot_state(model, states = c("YearlyAgeInc"), start_time = 59) + theme(legend.position = "none")